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ABSTRACT 


This  paper  describes  a  technique  for  the  solution  of  systems  of  symmetric  parabolic  partial 
differential  equations  that  are  associated  with  the  generation  of  biologic  textures  in  computer 
graphics.  The  technique  presented  allows  for  the  solution  of  the  texture  generating  equations  on 
the  highly  irregular  regions  associated  with  complex  image  geometries.  This  technique,  when 
used  with  an  automatic  grid  generation  program  allows  for  the  direct  and  efficient  generation  of 
textures  on  complex  surfaces.  Since  there  is  a  one  to  one  correspondence  between  the  generated 
texture  in  the  transformed  domain  and  the  original  surface  there  is  no  need  to  map  the  generated 
texture  pattern  back  to  the  surface  of  interest.  This  technique  can  be  applied  to  large  systems  of 
reaction  diffusion  equations,  thus  expanding  the  set  of  patterns  that  can  be  produced  using 
models  based  upon  reaction-diffusion  equations.  In  addition,  textures  based  upon  variable 
diffusion  constants  can  be  generated  including  a  wide  range  of  anisotropic  textures. 


KEYWORDS:  Texture  ;  Parabolic  Equations  ;  Reaction-Diffusion  ;  Irregular  Geometries  ; 
Iterative  methods ;  Mathematical  Modeling. 


OVERVIEW 

A  considerable  body  of  research  exists 
describing  the  generation  of  biologic  textures 
using  a  system  of  parabolic  partial  differential 
equations.  These  equations  are  called  reaction- 
diffusion  equations  in  applied  mathematics 
[Turk91]  [Witkin-Kass91].  Previous  work  in 
numerical  analysis  has  presented  techniques  to 
solve  such  equations  using  finite  difference 
techniques  on  problem  domains  that  are  square, 
circular  or  rectangular.  For  the  solution  of 
models  that  are  defined  upon  irregular  regions, 
previous  work  has  produced  techniques  that 
superimpose  square,  circular  or  rectangular 
grids  upon  the  irregular  problem  domain. 
Previous  work  has  solved  the  equations  on  these 


superimposed  grids  by  either  extrapolating 
boundary  values  to  grid  points  outside  the 
problem  domain  or  by  modifying  the  finite 
difference  stencils  used  in  approximating  the 
partial  derivatives  [Castillo  et  al.92]. 

This  paper  presents  a  technique  for  solving 
reaction-diffusion  equations  on  an  irregular 
region.  In  this  technique,  the  vertices  of  the 
image  are  transformed  from  physical  space  to  a 
new  coordinate  system  which  is  called  logical 
space.  The  system  of  reaction-diffusion 
equations  are  then  transformed  to  logical  space 
and  differenced  using  a  second  order  nearest 
neighbor  scheme.  Since  the  problem  domain  in 
logical  space  is  a  square  or  rectangle,  almost  any 
finite  difference  technique  developed  in  the  last 


reaction-diffusion  equations: 


100  years  can  be  used  to  produce  a  system  of 
difference  equations.  The  resulting  system  of 
difference  equations  is  solved  using  an  iterative 
scheme  to  produce  the  texture  values. 

The  primary  advantage  of  this  method  is  that  it 
eliminates  the  need  to  map  the  calculated  texture 
values  back  onto  the  original  untransformed 
surface  of  interest  Another  advantage  is  that  the 
method  is  very  easy  to  parallelize  for 
implementation  on  a  massively  parallel  computer 
architecture  such  as  a  Single  Instruction  Multiple 
Data  (SIMD)  architecture. 
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where  Lj  is  the  operator 


METHOD 

The  technique  described  in  this  paper  has  three 
goals  to  meet  for  the  generation  of  texture 
patterns:  grid  generation;  equation 
transformation;  and  solution  of  the  transformed 
system  of  equations. 

The  first  goal  of  this  method  is  to  generate  a 
grid.  To  accomplish  this  goal  the  coordinates  of 
the  vertices  of  the  image  to  be  textured  are 
transformed  from  physical  space  to  a 
rectangular  or  square  grid  in  logical  space.  This 
transformation  is  accomplished  using  Castillo's 
variational  grid  generation  method  [Castillo91] 
[Castillo-Richardson93].  This  grid  generation 
approach  controls  three  properties  of  the  grid: 
grid  spacing,  grid  cell  areas,  and  grid 
orthogonality.  The  grid  is  generated  by  direct 
minimization  of  a  discrete  functional.  By  grid 
orthogonality  is  meant  the  angles  of  the  lines 
between  the  lines  comprising  the  computational 
grid.  The  variational  grid  generation  method 
allows  the  user  to  vary  the  effects  of  above  three 
properties  upon  the  grid  generation  process. 
Grid  generation  provides  a  boundary  fitted 
coordinate  system  which  conforms  to  the 
geometry  of  the  region  where  the  PDE  or  system 
of  PDE's  are  specified. 

There  are  several  advantages  to  the  variational 
grid  generation  method.  The  weights  can  be 
chosen  to  prevent  folding  of  the  grid  and  to 
concentrate  grid  points  in  areas  of  the  grid  where 
there  are  large  variations  in  the  solution.  In 
addition,  the  generated  grid  can  be  adapted  to  the 
solution. 


f 


dX; 


)) 


The  rjj's  are  the  coefficients  of  the  reaction- 

diffusion  equation  and  consist  of  the  diffusion 
constants,  etc.  It  will  also  be  assumed  that  rjj  = 
rjj  so  that  Lj,.. ,  Ln  are  symmetric.  For  the  3-D 

case  the  i’s  and  j's  are  summed  from  1  to  3.  It 
should  be  noted  that  the  ry’s  can  be  constants  or 

variable.  The  system  of  equations  can  produce 
anisotropic  textures  subject  to  the  constraint  that 
the  operators  must  be  symmetric.  For  the 
situation  where  the  system  of  reaction-diffusion 
equations  are  linear  or  quasi  linear  this  is  not  a 
problem  since  these  types  of  systems  can  be 
modified  to  be  symmetric.  The  coefficients  of 
the  cross  derivatives  can  be  averaged  or 
similarly  modified  to  produce  symmetry  in  the 
system. 

Given  the  transformation  x  =  x(E)  where  E  =  (t, 
e  1 ,  e2  ,e3  )  we  have  the  following 

transformation  of  the  equations,  which  is  an 
extension  of  the  transformation  for  elliptic  partial 
differential  equations  given  by  [Steinberg91]  to 
the  case  of  parabolic  partial  differential  equations 
of  which  reaction-diffusion  equations  are  an 
example. 
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The  second  goal  of  this  method  is  to  transform 
the  system  of  reaction-diffusion  equations  to 
logical  space.  Consider  the  following  system  of 


ect  where  Aj  =  Eq.  (4) 
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In  order  to  preserve  the  original  properties  of  the 
system  of  reaction-diffusion  equations  the 
elements  of  the  cofactor  matrix  are  constrained 
to  satisfy  the  following  equation  called  the 
metric  identity.  The  metric  identity,  Eq.  (11), 
guarantees  the  second  order  approximation, 
symmetry  and  the  nearest  neighbor  property. 
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In  the  above  equations  J  is  the  Jacobian  of  the 
transformation  and  the  cjj's  are  the  entries  in  the 
cofactor  matrix  of  J  with  J  =  determinant(M). 
Note  that  the  sjj's  are  the  transformed 

coefficients. 

The  transformed  equations  can  then  be 
differenced  using  a  second  order  central 
difference  scheme.  Since  the  operators  of  the 
untransformed  equation  were  symmetric  that  the 
transformed  equations  are  symmetric  and  the 
difference  scheme  is  a  nearest  neighbor  scheme. 
It  should  be  noted  that  the  grid  generation 
method  provides  the  descretized  values  of  the 
transformed  coordinates.  For  example,  in  the 
two  dimensional  case,  we  have  the  following 
where  the  matrix  of  x  and  y  values  are  provided 
by  the  grid  generation  process.  Thus  the  partial 
derivatives  in  the  transformed  coordinate  system 
are  given  by  Eq.  (9)  and  Eq.  (10). 


If  the  grid  generation  method  transforms  the 
coordinates  of  the  comers  of  each  of  the  vertices 
in  rectangular  patches  comprising  the  image  then 
the  boundary  conditions  are  restricted.  If  the 
grid  generation  method  transforms  the 
coordinates  of  the  vertices  and  a  harmonic  or 
linear  average  of  the  comers  is  used  to  produce 
transformed  coordinates  at  the  transformed 
patch  centers  then  the  system  of  equations  can 
be  placed  in  conservative  form  and  a  much 
larger  class  of  boundary  conditions  can  be 
accommodated  including  neumann,  dirichlet  and 
robin  boundary  conditions.  The  Neumann  or 
no-flux  boundary  boundary  condition  is 
represented  by. 


(12)  Vn*/t  =  0 


The  Dirichlet  boundary  condition  is 
(13)  u  =  b(x,(E)) 

The  Robin  boundary  condition  is 


(9) 


_  (x(i  +  l,j)-x(i-lj)) 
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An  arbitrary  partial  differential  equation  is  said 
to  be  in  conservative  form  if  it  can  be  written  as 
the  divergence  of  a  vector  field.  Satisfaction  of 
the  metric  identity  implies  the  well  mannered 
behavior  of  the  conservative  form  and 
preservation  of  conservation  laws  throughout 
the  transformation  process.  The  vector  field  is 

called  the  flux  F  while  the  transformed  flux  is  d> 
and  they  arc  defined  by  the  following. 


(15) 
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The  system  of  transformed  and  differenced 
equations  are  then  solved  using  iterative 
methods  such  as  Successive  Over  Relaxation 
(SOR)  or  a  conjugate  gradient  iterative  scheme. 
Several  standard  textbooks  exist  on  the  subject 
of  SOR  and  conjugate  gradient  methods 
[Smith85]  [Strikwerda89],  It  should  be  noted 
that  the  conjugate  gradient  method  is  more 
amenable  to  parallel  computation  than  the  SOR 
method.  Preconditioning  of  the  conjugate 
gradient  scheme  is  recommended. 

The  metric  identity  guarantees  that  the  symmetry 
of  the  original  system  of  PDE's  is  preserved. 
Preserving  the  symmetry  of  the  original  system 
aids  in  producing  a  matrix  of  transformed 
coefficients  for  the  resulting  system  of 
difference  equations  that  is  well  banded  and  thus 
particularly  suited  for  efficient  implementation  of 
the  SOR  and  CG  iterative  methods. 

The  calculated  texture  values  are  the  same  in 
both  logical  and  physical  space.  This  is  a 
consequence  of  the  grid  generation  method  and 
the  metric  identity.  By  indexing  the  transformed 
points  in  logical  space  with  the  same  indices  as 
the  physical  space  the  need  to  map  the  texture 


values  back  onto  the  original  surface  is 
eliminated.  Thus,  the  warping  and 
discontinuities  produced  by  the  mapping  of 
textures  to  a  surface  by  algebraic  or  procedural 
techniques  are  avoided. 


TEXTURE  MODELING 

The  purpose  of  texture  generation  in  computer 
graphics  is  to  provide  texture  values  for  the 
geometrical  patches  that  comprise  an  image  in  2- 
D  or  3-D  space.  The  basic  assumption  is  that 
these  graphical  images  are  highly  irregular  in 
shape.  Because  of  the  irregular  nature  of  realistic 
graphics  images  a  large  repertoire  of  techniques 
has  been  developed  to  map  the  generated  texture 
values  onto  the  image  domain.  Several  examples 
are  parametric  texture  mapping  and  projective 
techniques  using  translations,  rotations  and 
standard  projections.  The  focus  of  this  paper  is 
to  provide  the  theoretical  framework  for  a  direct 
production  of  texture  values  for  a  computer 
image  free  of  the  necessity  to  map  the  generated 
texture  values  onto  the  image  geometry. 

Over  the  years  a  large  number  of  texture  models 
have  been  produced  using  algebraic  or 
procedural  techniques  to  generate  texture. 
Examples  are  marble  generation,  fractal  textures 
and  summation  of  sine  waves.  This  paper 
focuses  on  models  of  biologic  and  chemical 
processes  that  produce  patterns  suitable  for  use 
in  texture  synthesis. 

A  class  of  models  that  describe  pattern 
producing  processes  in  nature  consist  of 
systems  of  reaction-diffusion  equations. 
Reaction  diffusion  equations  are  parabolic  partial 
differential  equations  of  the  form  given  by  Eq. 
(1)  and  Eq.  (2)  above.  If  the  gi's  are  linear  then 
the  system  of  reaction-diffusion  equations  are 
linear.  If  the  gi's  are  nonlinear  then  the  system 
of  reaction-diffusion  equations  are  quasi  linear. 
Most  of  the  interesting  models  in  the  literature, 
such  as  Meinhardt’s  stripe  system  are  quasi 
linear. 

Reaction-diffusion  equations  are  examples  of 
phennomena  that  obey  conservation  laws.  If  we 
are  interested  in  a  concentration  or  texture  value 
at  a  point,  the  rate  of  change  of  this  texture 
value  is  related  to  the  inflow  and  outflow  of 
particles  (diffusion)  through  an  area  or  volume 
element  and  the  rate  of  particle  creation  or 
destruction  (reaction). 


There  are  many  examples  of  quasi  linear  models 
in  morphogenisis  and  population  biology  that 
embody  systems  of  reaetion-diffusion  equations 
and  that  form  patterns.  One  example  of  reaetion- 
diffusion  equations  that  form  beautiful  and 
interesting  patterns  is  the  model  of  reaction- 
diffusion  wave  propagation  based  upon  the 
eikonal  equation  which  has  the  following  form 
[Grindrod91].  This  model  produces  a  series  of 
spiral  and  toroidal  patterns  in  nature.  The  most 
interesting  patterns  are  produced  in  three 
dimensions.  Eq.  (17)  and  Eq.  (18)  describe 
these  spiral  and  toroidal  patterns. 


seconds  in  parallel  for  a  62  by  62  computational 
grid.  It  is  anticipated  that  the  timing  in  the 
parabolic  case  would  be  slightly  longer  for  each 
time  step. 

The  technique  presented  in  this  paper  is  a  direct 
method  for  solving  systems  of  symmetric 
parabolic  partial  differential  equations  on 
irregular  regions.  Application  of  this  technique 
to  systems  of  reaction-diffusion  equations  can 
be  used  to  produce  textures  in  a  straight  forward 
manner  free  from  the  problems  of  remapping 
back  onto  the  surface. 


(17)  A)+F(u)  =  0 

C/t  •  •  C/Xi  OXj 


The  wave  front  propagates  according  to  the 
following.  N  is  the  normal  to  the  wave  surface 
and  K  is  the  mean  curvature  of  the  wave. 
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CURRENT  AND  FUTURE  WORK 

Previous  work  using  this  technique  has  been 
performed  on  models  based  upon  elliptic  partial 
differential  equations.  Results  have  been 
obtained  for  the  transformation  of  elliptic  partial 
differential  equations  using  several  test  grids  and 
a  model  of  the  otolith  (eardrum)  membrane 
[Castillo  et  al.92]  [Tablewski90].  The  elliptic 
partial  differential  model  of  the  otolith  membrane 
has  also  been  implemented  in  parallel  on  a  mesh 
type  architecture  [Castillo-Richardson93]. 
Implementation  of  the  method  in  parallel  for  the 
case  of  elliptic  PDE's  resulted  in  a  reduction  of 
computation  time  by  a  factor  of  50-250. 

Current  work  involves  the  implementation  of  the 
grid  generation  method  in  parallel  on  mesh  type 
architecture  using  the  Applied  Memory 
Technology  DAP  610  computer.  In  addition, 
current  work  entails  the  implementation  of  the 
methods  presented  in  this  paper  to  the  case  of 
systems  of  parabolic  partial  differential 
equations. 
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